C $TITLE: 'RDPAT'
C $NOFLOATCALLS
C
      SUBROUTINE RDPAT(CUR,GAIN,AIR,AII,BIR,BII,CIR,CII,SI,CAB,SAB,
     1 BI,SALP,X,Y,Z,LD,LD3,LD4,IW)
C     COMPUTE RADIATION PATTERN, GAIN, NORMALIZED GAIN
CLARGE: CUR
      COMPLEX CUR
      COMPLEX*16 ETH,EPH,ERD
      COMPLEX*16 ZRATI,ZRATI2,T1,FRATI
      REAL*8 PI,TA,TD,IGNTP,IGAX,IGTP,HCIR,HBLK,HPOL,HCLIF,ISENS
      REAL*8 ATGN2,CANG,DB10,TSTOR1,TSTOR2,GCOP,GCON,ETHA,ETHM,ETHM2
      REAL*8 EPHA,EPHM,EPHM2,DFAZ,DFAZ2,CDFAZ,AIR,AII,BIR,BII,CIR,CII
      REAL*8 TILTA,STILTA,AXRAT,EMINR2,EMAJR2,GNMN,GNMJ
         REAL*8 GAIN
C     INTEGER HPOL,HBLK,HCIR,HCLIF
      INTEGER*4 N1,N2,N,NP,M1,M2,M,MP,IPSYM
      COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
      COMMON/SAVE/ KCOM,COM(19,5),EPSR,SIG,SCRWLT,SCRWRT,FMHZ
      COMMON/GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
     1 IFAR,IPERF,T1,T2
      COMMON/FPAT/NTH,NPH,IPD,IAVP,INOR,IAX,THETS,PHIS,DTH,DPH,RFLD,
     1 GNOR,CLT,CHT,EPSR2,SIG2,IXTYP,XPR6,PINR,PNLR,PLOSS,NEAR,NFEH,
     2 NRX,NRY,NRZ,XNR,YNR,ZNR,DXNR,DYNR,DZNR
      COMMON/PLOT/ IPLP1,IPLP2,IPLP3,IPLP4
      DIMENSION CUR(LD3),GAIN(LD4)
      DIMENSION CAB(LD),SAB(LD),SI(LD),BI(LD),SALP(LD),X(LD),Y(LD),Z(LD)
      DIMENSION AIR(LD),AII(LD),BIR(LD),BII(LD),CIR(LD),CII(LD)
      DIMENSION IGTP(4), IGAX(4), IGNTP(10), HPOL(3)
C***
      DATA HPOL/'LINEAR','RIGHT','LEFT'/,HBLK,HCIR/' ','CIRCLE'/
      DATA IGTP/'    - ','POWER ','- DIRE','CTIVE '/
      DATA IGAX/' MAJOR',' MINOR',' VERT.',' HOR. '/
      DATA IGNTP/' MAJOR',' AXIS ',' MINOR',' AXIS ','   VER','TICAL ','
     1 HORIZ','ONTAL ','      ','TOTAL '/
      DATA PI,TA,TD/3.141592654D0,1.745329252D-02,57.29577951D0/
C**
      DATA NORMAX/1200/
C**
      IF (IFAR.LT.2) GO TO 2
      WRITE(IW,35)
      IF (IFAR.LE.3) GO TO 1
      WRITE(IW,36)  NRADL,SCRWLT,SCRWRT
      IF (IFAR.EQ.4) GO TO 2
1     IF (IFAR.EQ.2.OR.IFAR.EQ.5) HCLIF=HPOL(1)
      IF (IFAR.EQ.3.OR.IFAR.EQ.6) HCLIF=HCIR
      CL=CLT/WLAM
      CH=CHT/WLAM
      ZRATI2=CDSQRT(1.D0/DCMPLX(EPSR2,-SIG2*WLAM*59.96))
      WRITE(IW,37)  HCLIF,CLT,CHT,EPSR2,SIG2
2     IF (IFAR.NE.1) GO TO 3
      WRITE(IW,41)
      GO TO 5
3     I=2*IPD+1
      J=I+1
      ITMP1=2*IAX+1
      ITMP2=ITMP1+1
      WRITE(IW,38)
      IF (RFLD.LT.1.E-20) GO TO 4
      EXRM=1./RFLD
      EXRA=RFLD/WLAM
      EXRA=-360.*(EXRA-AINT(EXRA))
      WRITE(IW,39)  RFLD,EXRM,EXRA
4     WRITE(IW,40)  IGTP(I),IGTP(J),IGAX(ITMP1),IGAX(ITMP2)
5     IF (IXTYP.EQ.0.OR.IXTYP.EQ.5) GO TO 7
      IF (IXTYP.EQ.4) GO TO 6
      PRAD=0.
      GCON=4.*PI/(1.+XPR6*XPR6)
      GCOP=GCON
      GO TO 8
6     PINR=394.51*XPR6*XPR6*WLAM*WLAM
7     GCOP=WLAM*WLAM*2.*PI/(376.73*PINR)
      PRAD=PINR-PLOSS-PNLR
      GCON=GCOP
      IF (IPD.NE.0) GCON=GCON*PINR/PRAD
8     I=0
      GMAX=-1.E10
      PINT=0.
      TMP1=DPH*TA
      TMP2=.5*DTH*TA
      PHI=PHIS-DPH
      DO 29 KPH=1,NPH
      PHI=PHI+DPH
      PHA=PHI*TA
      THET=THETS-DTH
      DO 29 KTH=1,NTH
      THET=THET+DTH
      IF (KSYMP.EQ.2.AND.THET.GT.90.01.AND.IFAR.NE.1) GO TO 29
      THA=THET*TA
      IF (IFAR.EQ.1) GO TO 9
      CALL FFLD (CUR,THA,PHA,ETH,EPH,X,Y,Z,SI,BI,
     1 SALP,AIR,AII,BIR,BII,CIR,CII,CAB,SAB,LD,LD3)
C**
      GO TO 10
9     CALL GFLD(RFLD/WLAM,PHA,THET/WLAM,ETH,EPH,ERD,ZRATI,KSYMP,LD,
     1 LD3,X,Y,Z,SI,BI,SALP,AIR,AII,BIR,BII,CIR,CII,CUR,CAB,SAB)
C**
C      ERDM=CABS(ERD)
      ERDM=ZABS(ERD)
      ERDA=CANG(ERD)
10    ETHM2=DREAL(ETH*DCONJG(ETH))
      ETHM=DSQRT(ETHM2)
      ETHA=CANG(ETH)
      EPHM2=DREAL(EPH*DCONJG(EPH))
      EPHM=DSQRT(EPHM2)
      EPHA=CANG(EPH)
      IF (IFAR.EQ.1) GO TO 28
C     ELLIPTICAL POLARIZATION CALC.
      IF (ETHM2.GT.1.E-20.OR.EPHM2.GT.1.E-20) GO TO 11
      TILTA=0.
      EMAJR2=0.
      EMINR2=0.
      AXRAT=0.
      ISENS=HBLK
      GO TO 16
11    DFAZ=EPHA-ETHA
      IF (EPHA.LT.0.) GO TO 12
      DFAZ2=DFAZ-360.
      GO TO 13
12    DFAZ2=DFAZ+360.
13    IF (ABS(DFAZ).GT.ABS(DFAZ2)) DFAZ=DFAZ2
      CDFAZ=DCOS(DFAZ*TA)
      TSTOR1=ETHM2-EPHM2
      TSTOR2=2.*EPHM*ETHM*CDFAZ
      TILTA=.5*ATGN2(TSTOR2,TSTOR1)
      STILTA=DSIN(TILTA)
      TSTOR1=TSTOR1*STILTA*STILTA
      TSTOR2=TSTOR2*STILTA*DCOS(TILTA)
      EMAJR2=-TSTOR1+TSTOR2+ETHM2
      EMINR2=TSTOR1-TSTOR2+EPHM2
      IF (EMINR2.LT.0.) EMINR2=0.
      AXRAT=DSQRT(EMINR2/EMAJR2)
      TILTA=TILTA*TD
      IF (AXRAT.GT.1.D-5) GO TO 14
      ISENS=HPOL(1)
      GO TO 16
14    IF (DFAZ.GT.0.) GO TO 15
      ISENS=HPOL(2)
      GO TO 16
15    ISENS=HPOL(3)
16    GNMJ=DB10(GCON*EMAJR2)
      GNMN=DB10(GCON*EMINR2)
      GNV=DB10(GCON*ETHM2)
      GNH=DB10(GCON*EPHM2)
      GTOT=DB10(GCON*(ETHM2+EPHM2))
      IF (INOR.LT.1) GO TO 23
      I=I+1
      IF (I.GT.NORMAX) GO TO 23
      GO TO (17,18,19,20,21), INOR
17    TSTOR1=GNMJ
      GO TO 22
18    TSTOR1=GNMN
      GO TO 22
19    TSTOR1=GNV
      GO TO 22
20    TSTOR1=GNH
      GO TO 22
21    TSTOR1=GTOT
22    GAIN(I)=TSTOR1
      IF (TSTOR1.GT.GMAX) GMAX=TSTOR1
23    IF (IAVP.EQ.0) GO TO 24
      TSTOR1=GCOP*(ETHM2+EPHM2)
      TMP3=THA-TMP2
      TMP4=THA+TMP2
      IF (KTH.EQ.1) TMP3=THA
      IF (KTH.EQ.NTH) TMP4=THA
C      DA=ABS(TMP1*(DCOS(TMP3)-DCOS(TMP4)))
      DA=ABS(TMP1*(COS(TMP3)-COS(TMP4)))
      IF (KPH.EQ.1.OR.KPH.EQ.NPH) DA=.5*DA
      PINT=PINT+TSTOR1*DA
      IF (IAVP.EQ.2) GO TO 29
24    IF (IAX.EQ.1) GO TO 25
      TMP5=GNMJ
      TMP6=GNMN
      GO TO 26
25    TMP5=GNV
      TMP6=GNH
26    ETHM=ETHM*WLAM
      EPHM=EPHM*WLAM
      IF (RFLD.LT.1.E-20) GO TO 27
      ETHM=ETHM*EXRM
      ETHA=ETHA+EXRA
      EPHM=EPHM*EXRM
      EPHA=EPHA+EXRA
27      CONTINUE
      WRITE(IW,42)  THET,PHI,TMP5,TMP6,GTOT,AXRAT,TILTA,ISENS,ETHM,ETHA
     1,EPHM,EPHA
C      GO TO 29
C***
C28    WRITE(IW,43)  RFLD,PHI,THET,ETHM,ETHA,EPHM,EPHA,ERDM,ERDA
      IF(IPLP1 .NE. 3) GO TO 299
      IF(IPLP3 .EQ. 0) GO TO 290
      IF(IPLP2 .EQ. 1 .AND. IPLP3 .EQ. 1) WRITE(8,*) THET,ETHM,ETHA
      IF(IPLP2 .EQ. 1 .AND. IPLP3 .EQ. 2) WRITE(8,*) THET,EPHM,EPHA
      IF(IPLP2 .EQ. 2 .AND. IPLP3 .EQ. 1) WRITE(8,*) PHI,ETHM,ETHA
      IF(IPLP2 .EQ. 2 .AND. IPLP3 .EQ. 2) WRITE(8,*) PHI,EPHM,EPHA
      IF(IPLP4 .EQ. 0) GO TO 299
290   IF(IPLP2 .EQ. 1 .AND. IPLP4 .EQ. 1) WRITE(8,*) THET,TMP5
      IF(IPLP2 .EQ. 1 .AND. IPLP4 .EQ. 2) WRITE(8,*) THET,TMP6
      IF(IPLP2 .EQ. 1 .AND. IPLP4 .EQ. 3) WRITE(8,*) THET,GTOT
C***
C***    ADD IPLP4.EQ.4 OPTION   RWA 03 APR 89   ADD 2 LINES
C***
      IF(IPLP2.EQ.1.AND.IPLP4.EQ.4) WRITE(8,*) THET,PHI,TMP5,TMP6,GTOT
      IF(IPLP2 .EQ. 2 .AND. IPLP4 .EQ. 1) WRITE(8,*) PHI,TMP5
      IF(IPLP2 .EQ. 2 .AND. IPLP4 .EQ. 2) WRITE(8,*) PHI,TMP6
      IF(IPLP2 .EQ. 2 .AND. IPLP4 .EQ. 3) WRITE(8,*) PHI,GTOT
C***
C***    ADD IPLP4.EQ.4 OPTION   RWA 03 APR 89   ADD 2 LINES
C***
      IF(IPLP2.EQ.1.AND.IPLP4.EQ.4) WRITE(8,*) PHI,THET,TMP5,TMP6,GTOT
      GO TO 299
28      CONTINUE
      WRITE(IW,43)  RFLD,PHI,THET,ETHM,ETHA,EPHM,EPHA,ERDM,ERDA
C***
C***    ADD IPLP2 = 3 FOR GROUND WAVE FIELDS  RWA 03 APR 89  ADD 11 LINES
C***
      IF(IPLP1.NE.3) GO TO 299
      IF (IPLP3.EQ.0) GO TO 299
      IF(IPLP2 .EQ. 1 .AND. IPLP3 .EQ. 1) WRITE(8,*) THET,ETHM,ETHA
      IF(IPLP2 .EQ. 1 .AND. IPLP3 .EQ. 2) WRITE(8,*) THET,EPHM,EPHA
      IF(IPLP2 .EQ. 1 .AND. IPLP3 .EQ. 3) WRITE(8,*) THET,ERDM,ERDA
      IF(IPLP2 .EQ. 2 .AND. IPLP3 .EQ. 1) WRITE(8,*) PHI,ETHM,ETHA
      IF(IPLP2 .EQ. 2 .AND. IPLP3 .EQ. 2) WRITE(8,*) PHI,EPHM,EPHA
      IF(IPLP2 .EQ. 2 .AND. IPLP3 .EQ. 3) WRITE(8,*) PHI,ERDM,ERDA
      IF(IPLP2 .EQ. 3 .AND. IPLP3 .EQ. 1) WRITE(8,*) RFLD,ETHM,ETHA
      IF(IPLP2 .EQ. 3 .AND. IPLP3 .EQ. 2) WRITE(8,*) RFLD,EPHM,EPHA
      IF(IPLP2 .EQ. 3 .AND. IPLP3 .EQ. 3) WRITE(8,*) RFLD,ERDM,ERDA
299   CONTINUE
C***
29    CONTINUE
      IF (IAVP.EQ.0) GO TO 30
      TMP3=THETS*TA
      TMP4=TMP3+DTH*TA*FLOAT(NTH-1)
C      TMP3=ABS(DPH*TA*FLOAT(NPH-1)*(DCOS(TMP3)-DCOS(TMP4)))
      TMP3=ABS(DPH*TA*FLOAT(NPH-1)*(COS(TMP3)-COS(TMP4)))
      PINT=PINT/TMP3
      TMP3=TMP3/PI
      WRITE(IW,44)  PINT,TMP3
30    IF (INOR.EQ.0) GO TO 34
      IF (ABS(GNOR).GT.1.E-20) GMAX=GNOR
      ITMP1=(INOR-1)*2+1
      ITMP2=ITMP1+1
      WRITE(IW,45)  IGNTP(ITMP1),IGNTP(ITMP2),GMAX
      ITMP2=NPH*NTH
      IF (ITMP2.GT.NORMAX) ITMP2=NORMAX
      ITMP1=(ITMP2+2)/3
      ITMP2=ITMP1*3-ITMP2
      ITMP3=ITMP1
      ITMP4=2*ITMP1
      IF (ITMP2.EQ.2) ITMP4=ITMP4-1
      DO 31 I=1,ITMP1
      ITMP3=ITMP3+1
      ITMP4=ITMP4+1
      J=(I-1)/NTH
      TMP1=THETS+FLOAT(I-J*NTH-1)*DTH
      TMP2=PHIS+FLOAT(J)*DPH
      J=(ITMP3-1)/NTH
      TMP3=THETS+FLOAT(ITMP3-J*NTH-1)*DTH
      TMP4=PHIS+FLOAT(J)*DPH
      J=(ITMP4-1)/NTH
      TMP5=THETS+FLOAT(ITMP4-J*NTH-1)*DTH
      TMP6=PHIS+FLOAT(J)*DPH
      TSTOR1=GAIN(I)-GMAX
      IF (I.EQ.ITMP1.AND.ITMP2.NE.0) GO TO 32
      TSTOR2=GAIN(ITMP3)-GMAX
      PINT=GAIN(ITMP4)-GMAX
31    WRITE(IW,46)  TMP1,TMP2,TSTOR1,TMP3,TMP4,TSTOR2,TMP5,TMP6,PINT
      GO TO 34
32    IF (ITMP2.EQ.2) GO TO 33
      TSTOR2=GAIN(ITMP3)-GMAX
      WRITE(IW,46)  TMP1,TMP2,TSTOR1,TMP3,TMP4,TSTOR2
      GO TO 34
33    WRITE(IW,46)  TMP1,TMP2,TSTOR1
34    RETURN
C
35    FORMAT (///,31X,39H- - - FAR FIELD GROUND PARAMETERS - - -,//)
36    FORMAT (40X,25HRADIAL WIRE GROUND SCREEN,/,40X,I5,6H WIRES,/,40X,1
     12HWIRE LENGTH=,F8.2,7H METERS,/,40X,12HWIRE RADIUS=,1P,E10.3,
     27H METERS)
37    FORMAT (40X,A6,6H CLIFF,/,40X,14HEDGE DISTANCE=,F9.2,7H METERS,/,4
     10X,7HHEIGHT=,F8.2,7H METERS,/,40X,15HSECOND MEDIUM -,/,40X,27HRELA
     2TIVE DIELECTRIC CONST.=,F7.3,/,40X,13HCONDUCTIVITY=,1P,E10.3,
     35H MHOS)
38    FORMAT(///,48X,30H- - - RADIATION PATTERNS - - -)
39    FORMAT (54X,6HRANGE=,1P,E13.6,7H METERS,/,54X,12HEXP(-JKR)/R=,
     1E12.5,9H AT PHASE,0P,F7.2,8H DEGREES,/)
40    FORMAT (/,2X,14H- - ANGLES - -,7X,2A6,7HGAINS -,7X,24H- - - POLARI
     1ZATION - - -,4X,20H- - - E(THETA) - - -,4X,16H- - - E(PHI) - -,2H
     2-,/,2X,5HTHETA,5X,3HPHI,7X,A6,2X,A6,3X,5HTOTAL,6X,5HAXIAL,5X,4HTIL
     3T,3X,5HSENSE,2(5X,9HMAGNITUDE,4X,6HPHASE ),/,2(1X,7HDEGREES,1X),3(
     46X,2HDB),8X,5HRATIO,5X,4HDEG.,8X,2(6X,7HVOLTS/M,4X,7HDEGREES))
41    FORMAT (///,28X,40H - - - RADIATED FIELDS NEAR GROUND - - -,//,8X,
     120H- - - LOCATION - - -,10X,16H- - E(THETA) - -,8X,14H- - E(PHI) -
     2 -,8X,17H- - E(RADIAL) - -,/,7X,3HRHO,6X,3HPHI,9X,1HZ,12X,3HMAG,6X
     3,5HPHASE,9X,3HMAG,6X,5HPHASE,9X,3HMAG,6X,5HPHASE,/,5X,6HMETERS,3X,
     47HDEGREES,4X,6HMETERS,8X,7HVOLTS/M,3X,7HDEGREES,6X,7HVOLTS/M,3X,7H
     5DEGREES,6X,7HVOLTS/M,3X,7HDEGREES,/)
42    FORMAT(1X,F7.2,F9.2,3X,3F8.2,F11.5,F9.2,2X,A6,2(1P,E15.5,0P,F9.2))
43    FORMAT (3X,F9.2,2X,F7.2,2X,F9.2,1X,3(3X,1P,E11.4,2X,0P,F7.2))
44    FORMAT (//,3X,19HAVERAGE POWER GAIN=,1P,E12.5,7X, 31HSOLID ANGLE U
     1SED IN AVERAGING=(,0P,F7.4,16H)*PI STERADIANS.,//)
45    FORMAT (//,37X,31H- - - - NORMALIZED GAIN - - - -,//,37X,2A6,4HGAI
     1N,/,38X,22HNORMALIZATION FACTOR =,F9.2,3H DB,//,3(4X,14H- - ANGLES
     2 - -,6X,4HGAIN,7X),/,3(4X,5HTHETA,5X,3HPHI,8X,2HDB,8X),/,3(3X,7HDE
     3GREES,2X,7HDEGREES,16X))
46    FORMAT (3(1X,2F9.2,1X,F9.2,6X))
      END
C
C
C
      SUBROUTINE FFLD(CUR,THET,PHI,ETH,EPH,X,Y,Z,SI,BI,
     1 SALP,AIR,AII,BIR,BII,CIR,CII,CAB,SAB,LD,LD3)
C
C     FFLD CALCULATES THE FAR ZONE RADIATED ELECTRIC FIELDS,
C     THE FACTOR EXP(J*K*R)/(R/LAMDA) NOT INCLUDED
C
      REAL*8 CONSX,ARG,DARG,PHX,PHY,ROX,ROY,ROZ,RRZ,ROZS,THX,THY,THZ,
     1 TTHET,D,DR,EL,A,B,C,OMEGA,SILL,TOP,BOT,BOO,TOO,RR,RI
      REAL*8 AIR,AII,BIR,BII,CIR,CII
CLARGE: CUR
      COMPLEX CUR
      COMPLEX*16 ZRATI,ZRATI2,T1,FRATI
      COMPLEX*16 CONST,EXA,GX,GY,GZ,CIX,CIY,CIZ,CCX,CCY,CCZ,CDP
      COMPLEX*16 ZRSIN,RRV,RRH,RRV1,RRH1,RRV2,RRH2,TIX,TIY,
     1 TIZ,ZSCRN,EX,EY,EZ,ETH,EPH
      INTEGER*4 N1,N2,N,NP,M1,M2,M,MP,IPSYM
      COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
      COMMON/GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
     1 IFAR,IPERF,T1,T2
      DIMENSION X(LD),Y(LD),Z(LD),SI(LD),BI(LD),SALP(LD)
      DIMENSION CAB(LD),SAB(LD),CUR(LD3),CONSX(2)
      DIMENSION AIR(LD),AII(LD),BIR(LD),BII(LD),CIR(LD),CII(LD)
      EQUIVALENCE (CONST,CONSX)
      DATA PI,TP,ETA/3.141592654D0,6.283185308D0,376.73/
      DATA CONSX/0.,-29.97922085D0/
C      PHX=-DSIN(PHI)
      PHX=-SIN(PHI)
C      PHY=DCOS(PHI)
      PHY=COS(PHI)
C      ROZ=DCOS(THET)
      ROZ=COS(THET)
      ROZS=ROZ
      THX=ROZ*PHY
      THY=-ROZ*PHX
C      THZ=-DSIN(THET)
      THZ=-SIN(THET)
      ROX=-THZ*PHY
      ROY=THZ*PHX
      IF (N.EQ.0) GO TO 20
C
C     LOOP FOR STRUCTURE IMAGE IF ANY
C
      DO 19 K=1,KSYMP
C
C     CALCULATION OF REFLECTION COEFFECIENTS
C
      IF (K.EQ.1) GO TO 4
      IF (IPERF.NE.1) GO TO 1
C
C     FOR PERFECT GROUND
C
      RRV=-(1.,0.)
      RRH=-(1.,0.)
      GO TO 2
C
C     FOR INFINITE PLANAR GROUND
C
1     ZRSIN=CDSQRT(1.-ZRATI*ZRATI*THZ*THZ)
      RRV=-(ROZ-ZRATI*ZRSIN)/(ROZ+ZRATI*ZRSIN)
      RRH=(ZRATI*ROZ-ZRSIN)/(ZRATI*ROZ+ZRSIN)
2     IF (IFAR.LE.1) GO TO 3
C
C     FOR THE CLIFF PROBLEM, TWO REFLCTION COEFFICIENTS CALCULATED
C
      RRV1=RRV
      RRH1=RRH
C      TTHET=DTAN(THET)
      TTHET=TAN(THET)
      IF (IFAR.EQ.4) GO TO 3
      ZRSIN=CDSQRT(1.-ZRATI2*ZRATI2*THZ*THZ)
      RRV2=-(ROZ-ZRATI2*ZRSIN)/(ROZ+ZRATI2*ZRSIN)
      RRH2=(ZRATI2*ROZ-ZRSIN)/(ZRATI2*ROZ+ZRSIN)
      DARG=-TP*2.*CH*ROZ
3     ROZ=-ROZ
      CCX=CIX
      CCY=CIY
      CCZ=CIZ
4     CIX=(0.,0.)
      CIY=(0.,0.)
      CIZ=(0.,0.)
C
C     LOOP OVER STRUCTURE SEGMENTS
C
      DO 17 I=1,N
      OMEGA=-(ROX*CAB(I)+ROY*SAB(I)+ROZ*SALP(I))
      EL=PI*SI(I)
      SILL=OMEGA*EL
      TOP=EL+SILL
      BOT=EL-SILL
      IF (ABS(OMEGA).LT.1.E-7) GO TO 5
      A=2.*DSIN(SILL)/OMEGA
      GO TO 6
5     A=(2.-OMEGA*OMEGA*EL*EL/3.)*EL
6     IF (ABS(TOP).LT.1.E-7) GO TO 7
      TOO=DSIN(TOP)/TOP
      GO TO 8
7     TOO=1.-TOP*TOP/6.
8     IF (ABS(BOT).LT.1.E-7) GO TO 9
      BOO=DSIN(BOT)/BOT
      GO TO 10
9     BOO=1.-BOT*BOT/6.
10    B=EL*(BOO-TOO)
      C=EL*(BOO+TOO)
      RR=A*AIR(I)+B*BII(I)+C*CIR(I)
      RI=A*AII(I)-B*BIR(I)+C*CII(I)
      ARG=TP*(X(I)*ROX+Y(I)*ROY+Z(I)*ROZ)
      IF (K.EQ.2.AND.IFAR.GE.2) GO TO 11
      EXA=DCMPLX(DCOS(ARG),DSIN(ARG))*DCMPLX(RR,RI)
C
C     SUMMATION FOR FAR FIELD INTEGRAL
C
      CIX=CIX+EXA*CAB(I)
      CIY=CIY+EXA*SAB(I)
      CIZ=CIZ+EXA*SALP(I)
      GO TO 17
C
C     CALCULATION OF IMAGE CONTRIBUTION IN CLIFF AND GROUND SCREEN
C     PROBLEMS.
C
11    DR=Z(I)*TTHET
C
C     SPECULAR POINT DISTANCE
C
      D=DR*PHY+X(I)
      IF (IFAR.EQ.2) GO TO 13
      D=DSQRT(D*D+(Y(I)-DR*PHX)**2)
      IF (IFAR.EQ.3) GO TO 13
      IF ((SCRWL-D).LT.0.) GO TO 12
C
C     RADIAL WIRE GROUND SCREEN REFLECTION COEFFICIENT
C
      D=D+T2
      ZSCRN=T1*D*DLOG(D/T2)
      ZSCRN=(ZSCRN*ZRATI)/(ETA*ZRATI+ZSCRN)
      ZRSIN=CDSQRT(1.-ZSCRN*ZSCRN*THZ*THZ)
      RRV=(ROZ+ZSCRN*ZRSIN)/(-ROZ+ZSCRN*ZRSIN)
      RRH=(ZSCRN*ROZ+ZRSIN)/(ZSCRN*ROZ-ZRSIN)
      GO TO 16
12    IF (IFAR.EQ.4) GO TO 14
      IF (IFAR.EQ.5) D=DR*PHY+X(I)
13    IF ((CL-D).LE.0.) GO TO 15
14    RRV=RRV1
      RRH=RRH1
      GO TO 16
15    RRV=RRV2
      RRH=RRH2
      ARG=ARG+DARG
16    EXA=DCMPLX(DCOS(ARG),DSIN(ARG))*DCMPLX(RR,RI)
C
C     CONTRIBUTION OF EACH IMAGE SEGMENT MODIFIED BY REFLECTION COEF. ,
C     FOR CLIFF AND GROUND SCREEN PROBLEMS
C
      TIX=EXA*CAB(I)
      TIY=EXA*SAB(I)
      TIZ=EXA*SALP(I)
      CDP=(TIX*PHX+TIY*PHY)*(RRH-RRV)
      CIX=CIX+TIX*RRV+CDP*PHX
      CIY=CIY+TIY*RRV+CDP*PHY
      CIZ=CIZ-TIZ*RRV
17    CONTINUE
      IF (K.EQ.1) GO TO 19
      IF (IFAR.GE.2) GO TO 18
C
C     CALCULATION OF CONTRIBUTION OF STRUCTURE IMAGE FOR INFINITE GROUND
C
      CDP=(CIX*PHX+CIY*PHY)*(RRH-RRV)
      CIX=CCX+CIX*RRV+CDP*PHX
      CIY=CCY+CIY*RRV+CDP*PHY
      CIZ=CCZ-CIZ*RRV
      GO TO 19
18    CIX=CIX+CCX
      CIY=CIY+CCY
      CIZ=CIZ+CCZ
19    CONTINUE
      IF (M.GT.0) GO TO 21
      ETH=(CIX*THX+CIY*THY+CIZ*THZ)*CONST
      EPH=(CIX*PHX+CIY*PHY)*CONST
      RETURN
20    CIX=(0.,0.)
      CIY=(0.,0.)
      CIZ=(0.,0.)
21    ROZ=ROZS
C
C     ELECTRIC FIELD COMPONENTS
C
      RFL=-1.
      DO 25 IP=1,KSYMP
      RFL=-RFL
      RRZ=ROZ*RFL
      CALL FFLDS(CUR(N+1),GX,GY,GZ,X,Y,Z,BI,ROX,ROY,RRZ,LD)
      IF (IP.EQ.2) GO TO 22
      EX=GX
      EY=GY
      EZ=GZ
      GO TO 25
22    IF (IPERF.NE.1) GO TO 23
      GX=-GX
      GY=-GY
      GZ=-GZ
      GO TO 24
23    RRV=CDSQRT(1.-ZRATI*ZRATI*THZ*THZ)
      RRH=ZRATI*ROZ
      RRH=(RRH-RRV)/(RRH+RRV)
      RRV=ZRATI*RRV
      RRV=-(ROZ-RRV)/(ROZ+RRV)
      ETH=(GX*PHX+GY*PHY)*(RRH-RRV)
      GX=GX*RRV+ETH*PHX
      GY=GY*RRV+ETH*PHY
      GZ=GZ*RRV
24    EX=EX+GX
      EY=EY+GY
      EZ=EZ-GZ
25    CONTINUE
      EX=EX+CIX*CONST
      EY=EY+CIY*CONST
      EZ=EZ+CIZ*CONST
      ETH=EX*THX+EY*THY+EZ*THZ
      EPH=EX*PHX+EY*PHY
      RETURN
      END
C
C
C
      SUBROUTINE GFLD (RHO,PHI,RZ,ETH,EPI,ERD,UX,KSYMP,LD,
     1 LD3,X,Y,Z,SI,BI,SALP,AIR,AII,BIR,BII,CIR,CII,CUR,CAB,SAB)
C
C     GFLD COMPUTES THE RADIATED FIELD INCLUDING GROUND WAVE.
C
      REAL*8 PI,TP,R,ARG,R1,R2,ZMH,ZPH,PHX,PHY,RX,RY,RIX,RIY,RIZ
      REAL*8 RHS,RHP,RHX,RHY,CALP,CBET,SBET,CPH,SPH,EL,RNX,RNY,RNZ
      REAL*8 RXYZ,OMEGA,SILL,TOP,BOT,BOO,TOO,A,B,C,RR,RI,THX,THY,THZ
      REAL*8 AIR,AII,BIR,BII,CIR,CII
CLARGE: CUR
      COMPLEX CUR
      COMPLEX*16 U,U2,XX1,XX2,EXA,ERV,EZV,ERH,EZH,EPH
      COMPLEX*16 CIX,CIY,CIZ,EX,EY,EPI,ETH,ERD,UX
      INTEGER*4 N1,N2,N,NP,M1,M2,M,MP,IPSYM
      COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
      COMMON/GWAV/ U,U2,XX1,XX2,R1,R2,ZMH,ZPH
      DIMENSION X(LD),Y(LD),Z(LD),SI(LD),BI(LD),SALP(LD)
      DIMENSION CAB(LD),SAB(LD),CUR(LD3)
      DIMENSION AIR(LD),AII(LD),BIR(LD),BII(LD),CIR(LD),CII(LD)
      DATA PI,TP/3.141592654D0,6.283185308D0/
C      R=DSQRT(RHO*RHO+RZ*RZ)
      R=SQRT(RHO*RHO+RZ*RZ)
      IF (KSYMP.EQ.1) GO TO 1
C      IF (CABS(UX).GT..5) GO TO 1
      IF (ZABS(UX).GT..5) GO TO 1
      IF (R.GT.1.D+5) GO TO 1
      GO TO 4
C
C     COMPUTATION OF SPACE WAVE ONLY
C
1     IF (RZ.LT.1.E-20) GO TO 2
C      THET=DATAN(RHO/RZ)
      THET=ATAN(RHO/RZ)
      GO TO 3
2     THET=PI*.5
3     CALL FFLD(CUR,THET,PHI,ETH,EPI,X,Y,Z,SI,BI,
     1 SALP,AIR,AII,BIR,BII,CIR,CII,CAB,SAB,LD,LD3)
      ARG=-TP*R
      EXA=DCMPLX(DCOS(ARG),DSIN(ARG))/R
      ETH=ETH*EXA
      EPI=EPI*EXA
      ERD=(0.,0.)
      RETURN
C
C     COMPUTATION OF SPACE AND GROUND WAVES.
C
4     U=UX
      U2=U*U
C      PHX=-DSIN(PHI)
      PHX=-SIN(PHI)
C      PHY=DCOS(PHI)
      PHY=COS(PHI)
      RX=RHO*PHY
      RY=-RHO*PHX
      CIX=(0.,0.)
      CIY=(0.,0.)
      CIZ=(0.,0.)
C
C     SUMMATION OF FIELD FROM INDIVIDUAL SEGMENTS
C
      DO 17 I=1,N
      DX=CAB(I)
      DY=SAB(I)
      DZ=SALP(I)
      RIX=RX-X(I)
      RIY=RY-Y(I)
      RHS=RIX*RIX+RIY*RIY
      RHP=DSQRT(RHS)
      IF (RHP.LT.1.E-6) GO TO 5
      RHX=RIX/RHP
      RHY=RIY/RHP
      GO TO 6
5     RHX=1.
      RHY=0.
6     CALP=1.-DZ*DZ
      IF (CALP.LT.1.E-6) GO TO 7
      CALP=DSQRT(CALP)
      CBET=DX/CALP
      SBET=DY/CALP
      CPH=RHX*CBET+RHY*SBET
      SPH=RHY*CBET-RHX*SBET
      GO TO 8
7     CPH=RHX
      SPH=RHY
8     EL=PI*SI(I)
      RFL=-1.
C
C     INTEGRATION OF (CURRENT)*(PHASE FACTOR) OVER SEGMENT AND IMAGE FOR
C     CONSTANT, SINE, AND COSINE CURRENT DISTRIBUTIONS
C
      DO 16 K=1,2
      RFL=-RFL
      RIZ=RZ-Z(I)*RFL
      RXYZ=DSQRT(RIX*RIX+RIY*RIY+RIZ*RIZ)
      RNX=RIX/RXYZ
      RNY=RIY/RXYZ
      RNZ=RIZ/RXYZ
      OMEGA=-(RNX*DX+RNY*DY+RNZ*DZ*RFL)
      SILL=OMEGA*EL
      TOP=EL+SILL
      BOT=EL-SILL
      IF(DABS(OMEGA).LT.1.D-7) GO TO 9
      A=2.*DSIN(SILL)/OMEGA
      GO TO 10
9     A=(2.-OMEGA*OMEGA*EL*EL/3.)*EL
10    IF(DABS(TOP).LT.1.D-7) GO TO 11
      TOO=DSIN(TOP)/TOP
      GO TO 12
11    TOO=1.-TOP*TOP/6.
12    IF(DABS(BOT).LT.1.D-7) GO TO 13
      BOO=DSIN(BOT)/BOT
      GO TO 14
13    BOO=1.-BOT*BOT/6.
14    B=EL*(BOO-TOO)
      C=EL*(BOO+TOO)
      RR=A*AIR(I)+B*BII(I)+C*CIR(I)
      RI=A*AII(I)-B*BIR(I)+C*CII(I)
      ARG=TP*(X(I)*RNX+Y(I)*RNY+Z(I)*RNZ*RFL)
      EXA=DCMPLX(DCOS(ARG),DSIN(ARG))*DCMPLX(RR,RI)/TP
      IF (K.EQ.2) GO TO 15
      XX1=EXA
      R1=RXYZ
      ZMH=RIZ
      GO TO 16
15    XX2=EXA
      R2=RXYZ
      ZPH=RIZ
16    CONTINUE
C
C     CALL SUBROUTINE TO COMPUTE THE FIELD OF SEGMENT INCLUDING
C     GROUND WAVE.
C
      CALL GWAVE (ERV,EZV,ERH,EZH,EPH)
      ERH=ERH*CPH*CALP+ERV*DZ
      EPH=EPH*SPH*CALP
      EZH=EZH*CPH*CALP+EZV*DZ
      EX=ERH*RHX-EPH*RHY
      EY=ERH*RHY+EPH*RHX
      CIX=CIX+EX
      CIY=CIY+EY
17    CIZ=CIZ+EZH
      ARG=-TP*R
      EXA=DCMPLX(DCOS(ARG),DSIN(ARG))
      CIX=CIX*EXA
      CIY=CIY*EXA
      CIZ=CIZ*EXA
      RNX=RX/R
      RNY=RY/R
      RNZ=RZ/R
      THX=RNZ*PHY
      THY=-RNZ*PHX
      THZ=-RHO/R
      ETH=CIX*THX+CIY*THY+CIZ*THZ
      EPI=CIX*PHX+CIY*PHY
      ERD=CIX*RNX+CIY*RNY+CIZ*RNZ
      RETURN
      END
C
C
C
      SUBROUTINE FFLDS(SCUR,EX,EY,EZ,XS,YS,ZS,S,ROX,ROY,ROZ,LD)
C     CALCULATES THE XYZ COMPONENTS OF THE ELECTRIC FIELD DUE TO
C     SURFACE CURRENTS
      REAL*8 TPI,CONSX,ARG,ROX,ROY,ROZ
CLARGE: SCUR
      COMPLEX SCUR
      COMPLEX*16 CONS,CT,EX,EY,EZ
      INTEGER*4 N1,N2,N,NP,M1,M2,M,MP,IPSYM
      COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
      DIMENSION XS(LD),YS(LD),ZS(LD),S(LD),CONSX(2),SCUR(1)
      EQUIVALENCE (CONS,CONSX)
      DATA TPI/6.283185308D0/,CONSX/0.,188.365D0/
      EX=(0.,0.)
      EY=(0.,0.)
      EZ=(0.,0.)
      I=LD+1
      DO 1 J=1,M
      I=I-1
      ARG=TPI*(ROX*XS(I)+ROY*YS(I)+ROZ*ZS(I))
      CT=DCMPLX(DCOS(ARG)*S(I),DSIN(ARG)*S(I))
      K=3*J
      EX=EX+SCUR(K-2)*CT
      EY=EY+SCUR(K-1)*CT
      EZ=EZ+SCUR(K)*CT
1     CONTINUE
      CT=ROX*EX+ROY*EY+ROZ*EZ
      EX=CONS*(CT*ROX-EX)
      EY=CONS*(CT*ROY-EY)
      EZ=CONS*(CT*ROZ-EZ)
      RETURN
      END
